Integrative Generalized Convex Clustering Optimization and Feature Selection for Mixed Multi-View Data

In mixed multi-view data, multiple sets of diverse features are measured on the same set of samples. By integrating all available data sources, we seek to discover common group structure among the samples that may be hidden in individualistic cluster analyses of a single data view. While several techniques for such integrative clustering have been explored, we propose and develop a convex formalization that enjoys strong empirical performance and inherits the mathematical properties of increasingly popular convex clustering methods. Specifically, our Integrative Generalized Convex Clustering Optimization (iGecco) method employs different convex distances, losses, or divergences for each of the different data views with a joint convex fusion penalty that leads to common groups. Additionally, integrating mixed multi-view data is often challenging when each data source is high-dimensional. To perform feature selection in such scenarios, we develop an adaptive shifted group-lasso penalty that selects features by shrinking them towards their loss-specific centers. Our so-called iGecco+ approach selects features from each data view that are best for determining the groups, often leading to improved integrative clustering. To solve our problem, we develop a new type of generalized multi-block ADMM algorithm using sub-problem approximations that more efficiently fits our model for big data sets. Through a series of numerical experiments and real data examples on text mining and genomics, we show that iGecco+ achieves superior empirical performance for high-dimensional mixed multi-view data.


Introduction
As the volume and complexity of data grows, statistical data integration has gained increasing attention as it can lead to discoveries which are not evident in analyses of a single data set. We study a specific data-integration problem where we seek to leverage common samples measured across multiple diverse sets of features that are of different types (e.g., continuous, count-valued, categorical, skewed continuous and etc.). This type of data is often called mixed, multi-view data (Hall and Llinas, 1997;Acar et al., 2011;Tang and Allen, 2018;Baker et al., 2019). While many techniques have been developed to analyze each individual data type separately, there are currently few methods that can directly analyze mixed multi-view data jointly. Yet, such data is common in many areas such as electronic health records, integrative genomics, multi-modal imaging, remote sensing, national security, online advertising, and environmental studies. For example in genomics, scientists often study gene regulation by exploring only gene expression data, but other data types, such as short RNA expression and DNA methylation, are all part of the same gene regulatory system. Joint analysis of such data can give scientists a more holistic view of the problem they study. But, this presents a major challenge as each individual data type is high-dimensional (i.e., a larger number of features than samples) with many uninformative features. Further, each data view can be of a different data type: expression of genes or short RNAs measured via sequencing is typically count-valued or zero-inflated plus skewed continuous data whereas DNA methylation data is typically proportion-valued. In this paper, we seek to leverage multiple sources of mixed data to better cluster the common samples as well as select relevant features that distinguish the inherent group structure.
We propose a convex formulation which integrates mixed types of data with different dataspecific losses, clusters common samples with a joint fusion penalty and selects informative features that separate groups. Due to the convex formulation, our methods enjoy strong mathematical and empirical properties. We make several methodological contributions. First, we consider employing different types of losses for better handling non-Gaussian data with Generalized Convex Clustering Optimization (Gecco), which replaces Euclidean distances in convex clustering with more general convex losses. We show that for different losses, Gecco's fusion penalty forms different types of centroids which we call loss-specific centers. To integrate mixed multi-view data and perform clustering, we incorporate different convex distances, losses, or divergences for each of the different data views with a joint convex fusion penalty that leads to common groups; this gives rise to Integrative Generalized Convex Clustering (iGecco). Further, when dealing with high-dimensional data, practitioners seek interpretability by identifying important features which can separate the groups. To facilitate feature selection in Gecco and iGecco, we develop an adaptive shifted group-lasso penalty that selects features by shrinking them towards their loss-specific centers, leading to Gecco+ and iGecco+ which performs clustering and variable selection simultaneously. To solve our problems in a computationally efficient manner, we develop a new general multiblock ADMM algorithm using sub-problem approximations, and make an optimization contribution by proving that this new class of algorithms converge to the global solution.

Related Literature
Our goal is to develop a unified, convex formulation of integrative clustering with feature selection based on increasingly popular convex clustering methods. Pelckmans et al. (2005); Lindsten et al. (2011);Hocking et al. (2011) proposed convex clustering which uses a fusion penalty to achieve agglomerative clustering like hierarchical clustering. This convex formulation guarantees a global optimal solution, enjoys strong statistical and mathematical theoretical properties, and often demonstrates superior empirical performance to competing approaches. Specifically, in literature, Pelckmans et al. (2005); Chi et al. (2017) showed it yields stable solutions to small perturbations on the data or tuning parameters; Radchenko and Mukherjee (2017) established clustering consistency by proving the clustering tree produced by convex clustering consistently estimates the clustering tree produced by the population procedure for the ℓ 1 penalty case; Tan and Witten (2015) established its link to hierarchical clustering as well as prediction consistency (finite sample bound for the prediction error); the perfect recovery properties of convex clustering with uniform weights have been proved by Zhu et al. (2014) for the two-clusters case and Panahi et al. (2017) for the general k-clusters case while Sun et al. (2021) proved results for general weighted convex clustering model; and many others have studied other appealing theoretical properties Chi and Steinerberger, 2019). Despite these advantages, convex clustering has not yet gained widespread popularity due to its intensive computation. Recently, some proposed fast and efficient algorithms to solve convex clustering and estimate its regularization paths (Chi and Lange, 2015;Weylandt et al., 2020). Meanwhile, convex clustering has been extended to biclustering (Chi et al., 2017) and many other applications Choi et al., 2019).
One potential drawback to convex clustering however, is that thus far, it has only been well-studied employing Euclidean distances between data points and their corresponding cluster centers. As is well known, the Euclidean metric suffers from poor performance with data that is highly non-Gaussian such as binary, count-valued, skewed data, or with data that has outliers. To alleviate the impact of outliers, Wang et al. (2016)  Despite these, however, no one has conducted a general investigation of convex clustering for non-Gaussian data, let alone studied data integration on mixed data, to the best of our knowledge. But, many others have proposed clustering methods for non-Gaussian data in other contexts. One approach is to perform standard clustering procedures on transformed data (Anders and Huber, 2010;Bullard et al., 2010;Marioni et al., 2008;Robinson and Oshlack, 2010). But, choosing an appropriate transformation that retains the original cluster signal is a challenging problem. Another popular approach is to use hierarchical clustering with specified distance metrics for non-Gaussian data (Choi et al., 2010;Fowlkes and Mallows, 1983). Closely related to this, Banerjee et al. (2005) studied different clustering algorithms utilizing a large class of loss functions via Bregman divergences. Yet, the proposed methods are all extensions of existing clustering approaches and hence inherit both good and bad properties of those approaches. There has also been work on model-based clustering, which assumes that data are generated by a finite mixture model; for example Banfield and Raftery (1993); Si et al. (2013) proposed such a model for the Poisson and negative binomial distributions. Still these methods have a non-convex formulation and local solutions like all model-based clustering methods. We propose to adopt the method similar to Banerjee et al. (2005) and study convex clustering using different loss functions; hence our method inherits the desirable properties of convex clustering and handles non-Gaussian data as well. More importantly, there is currently no literature on data integration using convex clustering and we achieve this by integrating different types of general convex losses with a joint fusion penalty.
Integrative clustering, however, has been well-studied in the literature. The most popular approach is to use latent variables to capture the inherent structure of multiple types of data. This achieves a joint dimension reduction and then clustering is performed on the joint latent variables (Shen et al., 2009Mo et al., 2013Mo et al., , 2017Meng et al., 2015). Similar in nature to the latent variables approach, matrix factorization methods assume that the data has an intrinsic low-dimensional representation, with the dimension often corresponding to the number of clusters Hellton and Thoresen, 2016;Zhang et al., 2012;Chalise and Fridley, 2017;Zhang et al., 2011;Yang and Michailidis, 2015). There are a few major drawbacks of latent variable or dimension reduction approaches, however. First it is often hard to directly interpret latent factors or low-dimensional projections. Second, many of these approaches are based on non-convex formulations yielding local solutions. And third, choosing the rank of factors or projections is known to be very challenging in practice and will often impact resulting clustering solutions. Another approach to integrative clustering is clustering of clusters (COC) which performs cluster analysis on every single data set and then integrates the primary clustering results into final group assignments using consensus clustering (Hoadley et al., 2014;Lock and Dunson, 2013;Kirk et al., 2012;Savage et al., 2013;Wang et al., 2014). This, however, has several potential limitations as each individual data set might not have enough signal to discern joint clusters or the individual cluster assignments are too disparate to reach a meaningful consensus. Finally, others have proposed to use distance-based clustering for mixed types of data by first defining an appropriate distance metric for mixed data (for example, the Gower distance by Gower, 1971) and then applying an existing distance-based clustering algorithm such as hierarchical clustering (Ahmad and Dey, 2007;Ji et al., 2012). Consequently, this method inherits both good and bad properties of distance-based clustering approaches. Notice that all of these approaches are either two-step approaches or are algorithmic or non-convex problems that yield local solutions. In practice, such approaches often lead to unreliable and unstable results.
penalized model-based clustering that selects features (Raftery and Dean, 2006;Wang and Zhu, 2008;Pan and Shen, 2007). Still, these methods inherit several advantages and disadvantages of model-based clustering approaches. Moreover, sparse integrative clustering is relatively under-studied. Shen et al. (2013); Mo et al. (2013) extended iCluster using a penalized latent variable approach to jointly model multiple omics data types. They induced sparsity on the latent variable coefficients via regularization. As feature selection is performed on the latent variables, however, this is less interpretable in terms of selecting features directly responsible for distinguishing clusters. Recently, and most closely related to our own work, Wang et al. (2018) proposed sparse convex clustering which adds a group-lasso penalty term on the cluster centers to shrink them towards zero, thus selecting relevant features. This penalty, however, is only appropriate for Euclidean distances when the data is centered; otherwise, the penalty term shrinks towards the incorrect cluster centers. For feature selection using different distances and losses, we propose an adaptive shifted group-lasso penalty that will select features by shrinking them towards their appropriate centroid.

Integrative Generalized Convex Clustering with Feature Selection
In this section, we introduce our new methods, beginning with the Gecco and iGecco and then show how to achieve feature selection via regularization. We also discuss some practical considerations for applying our methods and develop an adaptive version of our approaches.

Generalized Convex Clustering Optimization (Gecco)
In many applications, we seek to cluster data that is non-Gaussian. In the literature, most do this using different distance metrics other than Euclidean distances (Choi et al., 2010;Fowlkes and Mallows, 1983;de Souza and De Carvalho, 2004). Some use losses based on exponential family or deviances closely related to Bregman divergences (Banerjee et al., 2005).
To account for different types of losses for non-Gaussian data, we propose to replace the Euclidean distances in convex clustering with more general convex losses; this gives rise to Generalized Convex Clustering Optimization (Gecco): similar observations to be fused together simultaneously and is also rotation-invariant; but one could use ℓ 1 or ℓ ∞ -norm as well. γ is a positive tuning constant and w ij is a nonnegative weight. When γ equals zero, each data point occupies a unique cluster. As γ increases, the fusion penalty encourages some of the rows of the cluster center U to be exactly fused, forming clusters. When γ becomes sufficiently large, all centroids eventually coalesce to a single cluster centroid, which we define as the loss-specific center associated with ℓ ⋅ . Hence γ regulates both the cluster assignment and number of clusters, providing a family of clustering solutions. The weight w ij should be specified by the user in advance and is not a tuning parameter; we discuss choices of weights for various convex losses in Section 2.5.
Going beyond Euclidean distances, we propose to employ convex distance metrics as well as deviances associated with exponential family distributions and Bregman divergences, which are always convex. Interestingly, we show that each of these possible loss functions shrinks the cluster centers, U, to different loss-specific centers, instead of the mean-based centroid as in convex clustering with Euclidean distances. For example, one may want to use least absolute deviations (ℓ 1 -norm or Manhattan distances) for skewed data or for data with outliers; with this loss, we show that all observations fuse to the median when γ is sufficiently large. We emphasize loss-specific centers here as they will be important in feature selection in the next section. For completeness, we list common distances and deviance-based losses, as well as their loss-specific centers x j respectively in Table 1. (See Appendix F for all calculations associated with loss-specific centers, and we provide a formal proof when studying the properties of our approaches in Section 2.4.)

Integrative Generalized Convex Clustering (iGecco)
In data integration problems, we observe data from multiple sources and would like to get a holistic understanding of the problem by analyzing all the data simultaneously. In our framework, we integrate mixed multi-view data and perform clustering by employing different convex losses for each of the different data views with a joint convex fusion penalty that leads to common groups. Hence we propose Integrative Generalized Convex Clustering (iGecco) which can be formulated as follows: Here, we have K data sources. The k th data view X (k) is an n × p k matrix consisting of n observations and p k features; U (k) is also an n × p k matrix and the i th row, U i . k , is the cluster center associated with the point X i . k . And, ℓ k (X i . k , U i . k ) is the loss function associated with the k th data view. Still, we choose one loss type as ℓ k that is appropriate based on the data type of each view. Each loss function is weighted by π k , which is fixed by the user in advance. We have found that setting π k to be inversely proportional to the null deviance evaluated at the loss-specific center, i.e., π k = 1 ℓ k (X k , X k ) well in practice. The null deviance, ℓ k (X k , X k ) refers to the loss function evaluated at X k where each j th column of X k denotes the loss-specific center x j k . We employ this loss function weighting scheme to ensure equal scaling across data sets of different types. Recall that in generalized linear model (GLM), the likelihood-ratio test statistic, or the difference between the log-likelihoods, 2 ℓ k (X k , U k ) − ℓ k (X k , U 0 k ) , follows a χ 2 -distribution. Here U k is our iGecco estimate while U 0 k is the loss-specific center X k (cluster centroid when there is only one cluster). Therefore, the ratio of the two quantities, i.e., ℓ k (X k , U k ) ℓ k (X k , X k ) = π k ℓ k (X k , U k ) should be the same scale for each data view k. Finally, notice that we employ a joint convex fusion penalty on all of the U (k) 's; this incorporates information from each of the data sources and enforces the same group structure amongst the shared observations. Similar to convex clustering, our joint convex fusion penalty encourages the differences in the rows of concatenated centroids U 1 …U K to be shrunk towards zero, inducing a clustering behavior. Specifically, it forces the group structure of the i th row of U (k) to be the same for all k data views. Note the joint convex fusion penalty can be also written as . We say that subject i and i′ belong to the same cluster K , for all k. Hence, due to this joint convex fusion penalty, the common group structure property always holds. We study our methods further and prove some properties in Section 2.4.

Feature Selection: Gecco+ and iGecco+
In high dimensions, it is important to perform feature selection both for clustering purity and interpretability. Recently, Wang et al. (2018) proposed sparse convex clustering by imposing a group-lasso-type penalty on the cluster centers which achieves feature selection by shrinking noise features towards zero. This penalty, however, is only appropriate for Euclidean distances when the data is centered; otherwise, the penalty term shrinks towards the incorrect cluster centers. For example, the median is the cluster center with the ℓ 1 or Manhattan distances. Thus, to select features in this scenario, we need to shrink them towards the median, and we should enforce "sparsity" with respect to the median and not the origin. To address this, we propose adding a shifted group-lasso-type penalty which forces cluster center U ·j to shrink towards the appropriate loss-specific center x j for each feature.
Both the cluster fusion penalty and this new shifted-group-lasso-type feature selection penalty will shrink towards the same loss-specific center.
To facilitate feature selection with the adaptive shifted group-lasso penalty for one data type, our Generalized Convex Clustering Optimization with Feature Selection (Gecco+) is formulated as follows: Again, U is an n × p matrix and x j is the loss-specific center for the j th feature introduced in Table 1. The tuning parameter α controls the number of informative features and the feature weight ζ j is a user input which plays an important role to adaptively penalize the features.
(We discuss choices of ζ j in Section 2.5.2 when we introduce the adaptive version of our method.) When α is small, all features are selected and contribute to defining the cluster centers. When α grows sufficiently large, all features coalesce at the same value, the lossspecific center x j , and hence no features are selected and contribute towards determining the clusters. Another way of interpreting this is that the fusion penalty exactly fuses some of the rows of the cluster center U, hence determining groups of rows. On the other hand, the shifted group-lasso penalty shrinks whole columns of U towards their loss-specific centers, thereby essentially removing the effect of uninformative features. Selected features are then columns of U that are not shrunken to their loss-specific centers, U . j ≠ x j ⋅ 1 n .
These selected features, then, exhibit differences across the clusters determined by the fusion penalty. Clearly, sparse convex clustering of Wang et al. (2018) is a special case of Gecco+ using Euclidean distances with centered data. Our approach using both a row and column penalty is also reminiscent of convex biclustering (Chi et al., 2017) which uses fusion penalties on both the rows and columns to achieve checker-board-like biclusters.
Building upon integrative generalized convex clustering in Section 2.2 and our proposed feature selection penalty above, our Integrative Generalized Convex Clustering Optimization with Feature Selection (iGecco+) is formulated as follows: Again, U k is an n × p k matrix and x j k is the loss-specific center for the j th feature for k th data view. By construction, iGecco+ directly clusters mixed multi-view data and selects features from each data view simultaneously. Similarly, adaptive choice of ζ j k gives rise to adaptive iGecco+ which will be discussed in Section 2.5.2. Detailed discussions on practical choices of tuning parameters and weights can be also found in Section 2.5.

Properties
In this section, we develop some properties of our methods, highlighting several advantages of our convex formulation. Corresponding proofs can be found in Section A of the Appendix.
These properties illustrate some important advantages of our convex clustering approaches. Specifically, many other widely used clustering methods are known to suffer from poor local solutions, but any minimizer of our problem will achieve a global solution. Additionally, we show that iGecco+ is continuous with respect to the data, tuning parameters, and other input parameters. Together, these two properties are very important in practice and illustrate that the global solution of our method remains stable to small perturbations in the data and input parameters. Stability is a desirable property in practice as one would question the validity of a clustering result that can change dramatically with small changes to the data or parameters. Importantly, most popular clustering methods such as k-means, hierarchical clustering, model-based clustering, or low-rank based clustering, do not enjoy these same stability properties.
Finally in Proposition 3, we verify that when the tuning parameters are sufficiently large, full fusion of all observations to the loss-specific centers is achieved. Hence, our methods indeed behave as intended, achieving joint clustering of observations. We illustrate this property in Figure 1 where we apply Gecco+ to the authors data set (described fully in Section 4). Here, we illustrate how our solution, U γ, α , changes as a function of γ and α. This Similarly, when α is small, all features are selected and as α increases, some of the features get fused to their loss-specific center.

Practical Considerations and Adaptive iGecco+
In this section, we discuss some practical considerations for applying our method to real data. In the iGecco+ problem, π k , w and ζ j are user-specific fixed inputs while γ and α are tuning parameters; γ controls the number of clusters while α controls the number of features selected. We discuss choosing user-specific inputs such as weights as well as how to select tuning parameters. In doing so, we introduce an adaptive version of our method as well.

CHOICE OF WEIGHTS AND TUNING PARAMETERS-
In practice, a good choice of fusion weights w ij has been shown to enhance both computational efficiency and clustering quality of convex clustering (Chi and Lange, 2015). It has been empirically demonstrated that using weights inversely proportional to the distances yields superior clustering performance; this approach is widely adopted in practice. Further, setting many of the weights to zero helps reduce computation cost. Considering these two, the most common weights choice for convex clustering is to use K-nearest-neighbors method with a Gaussian kernel. Specifically, the weight between the sample pair (i, j) is set as w ij = I ij k exp −ϕd X i . , X j . , where I ij k equals 1 if observation j is among observation i's κ nearest neighbors or vice versa, and 0 otherwise. However, this choice of weights based on Euclidean distances may not work well for non-Gaussian data in Gecco(+) or for mixed data in iGecco(+). To account for different data types and better measure the similarity between observations, we still adopt K-nearest-neighbors method with an exponential kernel, but further extend this by employing appropriate distance metrics for specific data types in the exponential kernel. In particular, for weights in Gecco and Gecco+, we suggest using the same distance functions or deviances in the loss function of Gecco and Gecco+. For weights in iGecco and iGecco+, the challenge is that we need to employ a distance metric which measures mixed types of data. In this case, the Gower distance, which is a distance metric used to measure the dissimilarity of two observations measured in different data types (Gower, 1971), can address our problem. To be specific, the Gower distance between observation i and i′ overall is defined as d X i . , X i′ . = k 1 K j 1 refers to the Gower distance between observation i and i′ for feature j in data view k and R j k = max i, i′ |X ij k − X i′j k | is the range of feature j in data view k. In the literature, Gower distance has been commonly used as distance metrics for clustering mixed types of data (Wangchamhan et al., 2017;Hummel et al., 2017;Akay and Yüksel, 2018) and shown to yield superior performance than other distance metrics (Ali and Massmoudi, 2013;dos Santos and Zárate, 2015).
Alternatively, we also propose and explore using stochastic neighbor embedding weights based on symmetrized conditional probabilities (Maaten and Hinton, 2008). These have been shown to yield superior performance in high-dimensions and if there are potential outliers. Specifically, the symmetrized conditional probabilities are defined as p ij = p j i + p i j 2n , where p j i = exp(−ϕd(X i . , X j . )) k i exp( ϕd(X i X k )) . We propose to use the weights w ij = I ij k ⋅ p ij where I ij k still equals 1 if observation j is among observation i's κ nearest neighbors or vice versa, and 0 otherwise. Again, we suggest using distance metrics appropriate for specific data types or the Gower distance for mixed data. In empirical studies, we experimented with both weight choices and found that stochastic neighbor embedding weights tend to work better in high-dimensional settings and if there are outliers. Hence, we recommend these and employed them in our empirical investigations in Section 4 and 5.
Estimating the number of clusters in a data set is always challenging. Going beyond, we have two tuning parameters in our iGecco+ problem; γ controls the number of clusters while α controls the number of features selected. Current literature for tuning parameter selection for convex clustering mainly focuses on stability selection (Wang, 2010;Fang and Wang, 2012), hold-out validation (Chi et al., 2017) and information criterion (Tan and Witten, 2015). We first adopt the stability selection based approach for tuning parameter selection and follow closely the approach described in the work of Wang (2010); Fang and . We choose stability selection based approach because i) its selection consistency has been established and ii) Wang et al. (2018) adopted similar approach for tuning parameter selection and demonstrated strong empirical performance. However, stability selection is often computationally intensive in practice. To address this, we further explore information criterion based approaches like the Bayesian information criterion (BIC). We explain full details of both approaches in Appendix J and demonstrate empirical results when the number of clusters and features are not fixed but estimated based on the data.

ADAPTIVE GECCO+ AND IGECCO+ TO WEIGHT FEATURES-Finally,
we consider how to specify the feature weights, ζ j used in the shifted group-lasso penalty. While employing these weights are not strictly necessary, we have found, as did Wang et al. (2018), that like the fusion weights, well-specified ζ j 's can both improve performance and speed up computation. But unlike the fusion weights where we can use the pairwise distances, we don't have prior information on which features may be relevant in clustering. Thus, we propose to use an adaptive scheme that first fits the iGecco+ with no feature weights and uses this initial estimate to define feature importance for use in weights. This is similar to many adaptive approaches in the literature (Zou, 2006;Wang et al., 2018).
Our adaptive iGecco+ approach is given in Algorithm 1; this applies to adaptive Gecco+ as a special case as well. We assume that the number of clusters (or a range of the number of clusters) is known a priori. We begin by fitting iGecco+ with α = 1 and uniform feature weights ζ j k = 1. We then find the γ which gives the desired number of clusters, yielding the initial estimate, U k . (We provide alternative adaptive iGecco+ Algorithm 17 when the number of clusters is not known in Appendix J.) Next, similar to the adaptive approaches by Zou (2006); Wang et al. (2018), we use this initial estimate to adaptively weight features by proposing the following weights: ζ j ‖U . j k − x j k ⋅ 1 n ‖ 2 is close to zero in this case. Note, compared with sparse convex clustering where the authors defined feature weights ζ j by solving a convex clustering problem with feature penalty α = 0, we propose to fit iGecco+ with feature penalty α = 1 first and then update the feature weights adaptively. We find this weighting scheme works well in practice as it shrinks noise features more and hence penalizes more on those features. Such with-penalty initialization for adaptive weights has also been proposed in literature (Zhou et al., 2009;Fan et al., 2009;van de Geer et al., 2011). We also notice that noise features impact the distances used in the fusion weights as well. Hence, we suggest updating the distances adaptively by using the selected features to better measure the similarities between observations. To this end, we propose a new scheme to compute weighted Gower distances. First, we scale the features within each data view so that informative features in different data views contribute equally and on the same scale. Then, we employ the inverse of π k , i.e., the null deviance, to weight the distances from different data types, resulting in an aggregated and weighted Gower distance, d(X i . , X i′ . ) as further detailed in Algorithm 1.
Note that if the clustering signal from one particular data type is weak and there are few informative features for this data type, then our weighting scheme will down-weight this entire data type in the weighted Gower distance. In practice, our adaptive iGecco+ scheme works well as evidenced in our empirical investigations in the next sections.
2. Find γ which gives desired number of clusters; Get the estimate U k .
3. Update the feature weights ζ j k = 1/‖U . j k − x j k ⋅ 1 n ‖ 2 and fusion weights w ij = I ij κ exp(−ϕd(X i . , X i′ . )), where d(X i . , X i′ . ) = k 1 K j 1 p k ‖U j k x j k 1 n ‖ 2 max j ‖U j k x j k 1 n ‖ 2 ⋅ 1 π k ⋅ d ii′j k .
4. Fit adaptive iGecco+ with ζ , w and a sequence of γ and α; Find optimal γ and α which give desired number of clusters and features.

iGecco+ Algorithm
as a baseline approach. To save computation cost, we further develop a new multi-block ADMM-type procedure using inexact one-step approximation of the sub-problems. Our algorithm is novel from optimization perspective as we extend the multi-block ADMM to a higher number of blocks and combine it with the literature related to inexact-solve ADMM with sub-problem approximations, which often results in major computational savings.

Full ADMM to Solve iGecco+ (Naive Algorithm)
Given the shifted group-lasso and fusion penalties along with general losses, developing an optimization routine for iGecco+ method is less straight-forward than convex clustering or sparse convex clustering. In this section, we propose a simple ADMM algorithm to solve iGecco+ as a baseline algorithm and point out its drawbacks.
The most common approach to solve problems with more than two non-smooth functions is via multi-block ADMM (Lin et al., 2015;Deng et al., 2017), which decomposes the original problem into several smaller sub-problems and solves them in parallel at each iteration. Chen et al. (2016) established a sufficient condition for the convergence of three-block ADMM. We develop a multi-block ADMM approach to fit our problem for certain types of losses and prove its convergence.
We first recast iGecco+ problem (1) as the equivalent constrained optimization problem: Recently, Weylandt et al. (2020) derived the ADMM for convex clustering in matrix form and we adopt similar approach. We index a centroid pair by l = l 1 , l 2 with l 1 < l 2 , define the set of edges over the non-zero weights ε = l = l 1 , l 2 : w l > 0 , and introduce a new variable k to account for the difference between the two centroids. Hence V k is a matrix containing the pairwise differences between connected rows of U k and the constraint is equivalent to DU k − V k = 0 for all k; D ∈ ℝ ε × n is the directed difference matrix corresponding to the non-zero fusion weights. We give general-form multi-block ADMM (Algorithm 2) to solve iGecco+. Here prox ℎ ⋅ x = argmin z 1 2 x − z 2 2 + ℎ z is the proximal mapping of function h. Also, the superscript in U k in Algorithm 2 refers to the k th data view; we omit iteration counter indices in all iGecco+ algorithm for notation purposes and use the most recent values of the updates. The dual variable is denoted by Λ k .

Algorithm 2
General multi-block algorithm for iGecco+ Notice that, in Algorithm 2, there is no closed-form analytical solution for the U (k) subproblem for general losses. Typically, at each iteration of Algorithm 2, we need to apply an inner optimization routine, which requires a nested iterative solver, to solve the U (k) sub-problem until full convergence. In the next section, we seek to speed up this algorithm by using U (k) sub-problem approximations. But, first we propose two different approaches to fully solve the U (k) sub-problem based on specific loss types and then use these to develop a one-step update to solve the sub-problem approximately with guaranteed convergence. For the V sub-problem, one can easily show that it has a closed-form analytical solution for each iteration, as given in Algorithm 2.

iGecco+ Algorithm
We have introduced Algorithm 2, a simple baseline ADMM approach to solve iGecco+. In this section, we consider different ways to solve the U (k) sub-problem in Algorithm 2. First, based on specific loss types (differentiable or non-differentiable), we propose two different algorithms to solve the U (k) sub-problem to full convergence. These approaches, however, are rather slow for general losses as there is no closed-form solution which requires another nested iterative solver. To address this and in connection with current literature on variants of ADMM with sub-problem approximations, we propose iGecco+ algorithm, a multi-block ADMM algorithm which solves the sub-problems approximately by taking a single one-step update. We prove convergence of this general class of algorithms, a novel result in the optimization literature.

DIFFERENTIABLE CASE-When
the loss ℓ k is differentiable, we consider solving the U (k) sub-problem with proximal gradient descent, which is often used when the objective function can be decomposed into a differentiable and a non-differentiable function. While there are many other possible optimization routines to solve the U (k) sub-problem, we choose proximal gradient descent as there is existing literature proving convergence of ADMM algorithms with approximately solved sub-problems using proximal gradient descent Lu et al., 2016). We will discuss in detail how to approximately solve the sub-problem by taking a one-step approximation in Section 3.2.3. Based upon this, we propose Algorithm 3, which solves the U (k) sub-problem by running full iterative proximal gradient descent to convergence. Here P 2 (U k ; ζ k ) = j 1 p k ζ j k ‖U j k ‖ 2 .

Algorithm 3
U (k) sub-problem for differentiable loss ℓ k (proximal gradient): In Algorithm 3 and typically in general (proximal gradient) descent algorithms, we need to choose an appropriate step size s k to ensure convergence. Usually we employ a fixed step size by computing the Lipschitz constant as in the squared error loss case; but in our method, it is hard to compute the Lipschitz constant for most of our general losses. Instead, we suggest using backtracking line search procedure proposed by Beck and Teboulle (2009); Parikh et al. (2014), which is a common way to determine step size with guaranteed convergence in optimization. Further, we find decomposing the U (k) sub-problem to p k separate U . j k sub-problems brings several advantages such as i) better convergence property (than updating U (k) 's all together) due to adaptive step size for each U . j k sub-problem and ii) less computation cost by solving each in parallel. Hence, in this case, we propose to use proximal gradient for each separate U . j k sub-problem. To achieve this, we assume that the loss is elementwise, which is satisfied by every deviance-based loss. Last, as mentioned, there are many other possible ways to solve the U (k) sub-problem than proximal gradient, such as ADMM. We find that when the loss function is squared Euclidean distances or the loss function has a Hessian matrix that can be upper bounded by a fixed matrix, using ADMM approach saves more computation. We provide all implementation details discussed above in Section C of the Appendix.

NON-DIFFERENTIABLE CASE-When
the loss ℓ k is non-differentiable, we can no longer adopt the proximal gradient method to solve the U (k) sub-problem as the objective is now a sum of more than one separable non-smooth function. To address this, as mentioned, we can use multi-block ADMM; in this case, we introduce new blocks for the non-smooth functions and hence develop a full three-block ADMM approach to fit our problem.
To augment the non-differentiable term, we assume that our loss function can be written as ℓ k (X (k) , U (k) ) = f k (g k (X (k) , U (k) )) where f k is convex but non-differentiable and g k is affine.
This condition is satisfied by all distance-based losses with g k (X (k) , U (k) ) = X (k) − U (k) ; for example, for Manhattan distances, we have f k (Z) = j 1 p z j 1 vec(Z) 1 , and g k (X, U) = X − U. The benefit of doing this is that now the U (k) sub-problem has a closed- where X k is an n × p k matrix with the j th column equal to scalar x j k .
It is clear that we can use multi-block ADMM to solve the problem above and each primal variable has a simple update with a closed-form solution. We propose Algorithm 4, a full, iterative multi-block ADMM, to solve the U (k) sub-problem when the loss is a non-differentiable distance-based function. Algorithm 4 applies to iGecco+ with various distances such as Manhattan, Minkowski and Chebychev distances and details are given in Section D of the Appendix.

IGECCO+ ALGORITHM: FAST ADMM WITH INEXACT ONE-STEP APPROXIMATION TO THE SUB-
PROBLEM-Notice that for both Algorithm 3 and 4, we need to run them iteratively to full convergence in order to solve the U (k) sub-problem in Algorithm 2 for each iteration, which is dramatically slow in practice. To address this, in literature, many have proposed variants of ADMM with guaranteed convergence that find an inexact, one-step, approximate solution to the sub-problem (without fully solving it); these include the generalized ADMM (Deng and Yin, 2016), proximal ADMM (Shefi and Teboulle, 2014; Banert et al., 2016) and proximal linearized ADMM Lu et al., 2016). Thus, we propose to solve the U (k) sub-problem approximately by taking a single one-step update of the algorithm (Algorithm 3 or 4) for both types of losses and prove convergence. For the differentiable loss case, we propose to apply the proximal linearized ADMM approach while for the non-differentiable case, we show that taking a one-step update of Algorithm 4, along with V and Λ update in Algorithm 2, is equivalent to applying a four-block ADMM to the original iGecco+ problem and we provide a sufficient condition for the convergence of four-block ADMM. Our algorithm, to the best of our knowledge, is the first to incorporate higher-order multi-block ADMM and inexact ADMM with a one-step update to solve sub-problems for general loss functions.
When the loss is differentiable, as mentioned in Algorithm 3, one can use full iterative proximal gradient to solve the U . j k sub-problem, which however, is computationally burdensome. To avoid this, many proposed variants of ADMM which find approximate solutions to the sub-problems. Specifically, closely related to our problem here, Liu et al. (2013); Lu et al. (2016) proposed proximal linearized ADMM which solves the subproblems efficiently by linearizing the differentiable part and then applying proximal gradient due to the non-differentiable part. We find their approach fits into our problem and hence develop a proximal linearized 2-block ADMM to solve iGecco+ when the loss ℓ k is differentiable and gradient is Lipschitz continuous. It can be shown that applying proximal linearized 2-block ADMM to Algorithm 2 is equivalent to taking a one-step update of Algorithm 3 along with V and Λ update in Algorithm 2. In this way, we avoid running full iterative proximal gradient updates to convergence for the U (k) sub-problem as in Algorithm 3 and hence save computation cost.
When the loss is non-differentiable, we still seek to take an one-step update to solve the U (k) sub-problem. We achieve this by noticing that taking a one-step update of Algorithm 4 along with V and Λ update in Algorithm 2 is equivalent to applying multi-block ADMM to the original iGecco+ problem recast as follows (for non-differentiable distance-based loss): Typically, general higher-order multi-block ADMM algorithms do not always converge, even for convex functions (Chen et al., 2016). We prove convergence of our algorithm and establish a novel convergence result by casting the iGecco+ with non-differentiable losses as a four-block ADMM, proposing a sufficient condition for convergence of higher-order multi-block ADMMs, and finally showing that our problem satisfies this condition. (Details are given in the proof of Theorem 4 in Appendix B.) Therefore, taking a one-step update of Algorithm 4 converges for iGecco+ with non-differentiable losses.
So far, we have proposed inexact-solve one-step update approach for both differentiable loss and non-differentiable loss case. For mixed type of losses, we combine these two algorithms and this gives Algorithm 5, a multi-block ADMM algorithm with inexact one-step approximation to the U (k) sub-problem to solve iGecco+. We also establish the following convergence result.
Theorem 4 (iGecco+ convergence) Consider the iGecco+ problem (1) with fixed inputs π k , w and ζ j . If ℓ k is convex for all k, Algorithm 5 converges to a global solution. In addition, if each ℓ k is strictly convex, it converges to the unique global solution.
Remark. Our corresponding Theorem 4 establishes a novel convergence result as it is the first to show the convergence of four-block or higher ADMM using approximate subproblems for both differentiable and non-differentiable losses.
It is easy to see that Algorithm 5 can be applied to solve other Gecco-related methods as special cases. When K = 1, Algorithm 5 gives the algorithm to solve Gecco+. When α = 0, Algorithm 5 gives the algorithm to solve iGecco+. When K = 1 and α = 0, Algorithm 5 gives the algorithm to solve Gecco.
To conclude this section, we compare the convergence results of using both full ADMM and inexact ADMM with one-step update in the sub-problem to solve Gecco+ (n = 120 and p = 210) in Figure 2. The left plots show the number of iterations needed to yield optimization convergence while the right plots show computation time. We see that Algorithm 5 (one-step update to solve the sub-problem) saves much more computational time than Algorithm 2 (full updates to solve the sub-problem). It should be pointed out that though Algorithm 5 takes more iterations to converge due to inexact approximation for each iteration, we still reduce computation time dramatically as the computation time per iteration is much less than the full-solve approach.

Simulation Studies
In this section, we first evaluate the performance of Gecco+ against existing methods on non-Gaussian data. Next we compare iGecco+ with other methods on mixed multi-view data.

Non-Gaussian Data
In this subsection, we evaluate the performance of Gecco and (adaptive) Gecco+ by comparing it with k-means, hierarchical clustering and sparse convex clustering. For simplicity, we have the following naming convention for all methods: loss type name + Gecco(+). For example, Poisson Deviance Gecco+ refers to Generalized Convex Clustering with Feature Selection using Poisson deviance. Sparse CC refers to sparse convex clustering using Euclidean distances where each feature is centered first. We measure the accuracy of clustering results using adjusted Rand index (Hubert and Arabie, 1985). The adjusted Rand index is the corrected-for-chance version of the Rand index, which is used to measure the agreement between the estimated clustering assignment and the true group label. A larger adjusted Rand index implies a better clustering result. For all methods we consider, we assume oracle number of clusters for fair comparisons.
Each simulated data set is comprised of n = 120 observations with 3 clusters. Each cluster has an equal number of observations. Only the first 10 features are informative while the rest are noise. We consider the following simulation scenarios.
• S1: Spherical data with outliers The first 10 informative features in each group are generated from a Gaussian distribution with different μ k 's for each class. Specifically, the first 10 features are generated from N(μ k , I 10 ) where μ 1 = ( − 2.5 · 1 5 T . The outliers in each class are generated from a Gaussian distribution with the same mean centroid μ k but with higher variance, i.e., N(μ k , 25 · I 10 ). The remaining noise features are generated from N(0,1).
In the first setting (S1A), number of noise features ranges in 25, 50, 75, ⋯ up to 225 with the proportion of number of outliers fixed ( = 5%). We also consider the setting when the variance of noise features increases with number of features fixed p = 200 and number of outliers fixed (S1B) and high-dimensional setting where p ranges from 250, 500, 750 to 1000 (S1C).
• S2: Non-spherical data with three half moons Here we consider the standard simulated data of three interlocking half moons as suggested by Chi and Lange (2015) and Wang et al. (2018). The first 10 features are informative in which each pair makes up two-dimensional three interlocking half moons. We randomly select 5% of the observations in each group and make them outliers. The remaining noise features are generated from N(0,1). The number of noise features ranges from 25, 50, 75, ⋯ up to 225. In both S1 and S2, we compare Manhattan Gecco+ with other existing methods.
We summarize simulation results in Figure 3. We find that for spherical data with outliers, adaptive Manhattan Gecco+ performs the best in high dimensions. Manhattan Gecco performs well in low dimensions but poorly as the number of noisy features increases. Manhattan Gecco+ performs well as the dimension increases, but adaptive Manhattan Gecco+ outperforms the former as it adaptively penalizes the features, meaning that noisy features quickly get zeroed out in the clustering path and that only the informative features perform important roles in clustering. We see that, without adaptive methods, we do not achieve the full benefit of performing feature selection. As we perform adaptive Gecco+, we show vast improvement in clustering purity as the number of noise features grows where regular Gecco performs poorly. Sparse convex clustering performs the worst as it tends to pick outliers as singleton clusters. In the presence of outliers, Manhattan Gecco+ performs much better than sparse convex clustering as we choose a loss function that is more robust to outliers. Interestingly, k-means performs better than sparse convex clustering. This is mainly because sparse convex clustering calculates pairwise distance in the weights, placing outliers in singleton clusters more likely than k-means which calculates within-cluster variances (where outliers could be closer to the cluster centroids than to other data points). Our simulation results also show that adaptive Manhattan Gecco+ works well for non-spherical data by selecting the correct features. For count data, all three adaptive Gecco+ methods perform better than k-means, hierarchical clustering and sparse convex clustering. We should point out that there are several linkage options for hierarchical clustering. For visualization purposes, we only show the linkage with the best and worst performance instead of all the linkages. Also we use the appropriate data-specific distance metrics in hierarchical clustering. For k-means, we use k-means++ (Arthur and Vassilvitskii, 2006) for initialization. Table 2 shows the variable selection accuracy of sparse convex clustering and adaptive Gecco+ in terms of F 1 score. In all scenarios, we fix p = 225. We see that adaptive Gecco+ selects the correct features, whereas sparse convex clustering performs poorly.

Multi-View Data
In this subsection, we evaluate the performance of iGecco and (adaptive) iGecco+ on mixed multi-view data by comparing it with hierarchical clustering, iClusterPlus  and Bayesian Consensus Clustering (Lock and Dunson, 2013). Again, we measure the accuracy of clustering results using the adjusted Rand index (Hubert and Arabie, 1985).
As before, each simulated data set is comprised of n = 120 observations with 3 clusters.
Each cluster has an equal number of observations. Only the first 10 features are informative while the rest are noise. We have three data views consisting of continuous data, countvalued data and binary/proportion-valued data. We investigate different scenarios with different dimensions for each data view and consider the following simulation scenarios: • S1: Spherical data with p 1 = p 2 = p 3 = 10 • S2: Three half-moon data with p 1 = p 2 = p 3 = 10 We employ a similar simulation setup as in Section 4.1 to generate each data view. The difference is that here for informative features, we increase the within-cluster variance for Gaussian data and decrease difference of cluster mean centroids μ k 's for binary and count data so that there is overlap between different clusters. Specifically, for spherical cases, Gaussian data is generated from N(μ k , 3 · I 10 ); count data is generated from Poisson with different μ k 's (μ 1 = 2, μ 2 = 4, μ 3 = 6, etc); binary data is generated from Bernoulli with different μ k 's (μ 1 = 0.5, μ 2 = 0.2, μ 3 = 0.8, etc). For half-moon cases, continuous data is simulated with larger noise and the count and proportion-valued data is generated via a copula transform. In this manner, we have created a challenging simulation scenario where accurate clustering results cannot be achieved by considering only a single data view.
Again, for fair comparisons across methods, we assume oracle number of clusters. When applying iGecco(+) methods, we employ Euclidean distances for continuous data, Manhattan distances for count-valued data and Bernoulli log-likelihood for binary or proportion-valued data. We use the latter two losses as they perform well compared with counterpart losses in Gecco+ and demonstrate faster computation speed.
Simulation results in Table 3 and Table 4 show that our methods perform better than existing methods. In low dimensions, iGecco performs comparably with iCluster and Bayesian Consensus Clustering for spherical data. For non-spherical data, iGecco performs much better. For high dimensions, iGecco+ performs better than iGecco while adaptive iGecco+ performs the best as it achieves the full benefit of feature selection. We also applied k-means to the simulated data. The results of k-means are close to (or in some cases worse than) hierarchical clustering with the best-performing linkage; hence we only show the results of hierarchical clustering here for comparison.
Also we show the variable selection results in Table 5 and compare our method to that of iClusterPlus. For fair comparisons, we assume oracle number of features. For our method, we choose α that gives oracle number of features; for iClusterPlus, we select top features based on Lasso coefficient estimates. Our adaptive iGecco+ outperforms iClusterPlus for all scenarios.
Note in this section, we assume that the oracle number of clusters and features are known a priori for fair comparisons. Results when the number of clusters and features are not fixed but estimated based on the data using tuning parameter selection, are given in Appendix J.3.

Real Data Examples
In this section, we apply our methods to various real data sets and evaluate our methods against existing ones. We first evaluate the performance of Gecco+ for several real data sets and investigate the features selected by various Gecco+ methods.

Authors Data
The authors data set consists of word counts from n = 841 chapters written by four famous English-language authors (Austen, London, Shakespeare, and Milton). Each class contains an unbalanced number of observations with 69 features. The features are common "stop words" like "a", "be" and "the" which are typically removed before text mining analysis. We use Gecco+ not only to cluster book chapters and compare the clustering assignments with true labels of authors, but also to identify which key words help distinguish the authors. We choose tuning parameters using BIC based approach; results for stability selection based approach are given in Table 15, Appendix J.3.
In Table 6, we compare Gecco+ with existing methods in terms of clustering quality. For hierarchical clustering, we only show the linkage with the best performance (in this whole section). Our method outperforms k-means and the best hierarchical clustering method. Secondly, we look at the word texts selected by Gecco+. As shown in Table 7, Jane Austen tended to use the word "her" more frequently than the other authors; this is expected as the subjects of her novels are typically females. The word "was" seems to separate Shakespeare and Jack London well. Shakespeare preferred to use present tense more while Jack London preferred to use past tense more. To summarize, our Gecco+ not only has superior clustering performance but also selects interpretable features.

TCGA Breast Cancer Data
The From Table 8, our method outperforms k-means and the best hierarchical clustering method. Next, we look at the genes selected by Gecco+ in Table 9. FOXA1 is known to be a key gene that characterizes luminal subtypes in DNA microarray analyses (Badve et al., 2007). GATA binding protein 3 (GATA3) is a transcriptional activator highly expressed by the luminal epithelial cells in the breast (Mehra et al., 2005). ERBB2 is known to be associated with HER2 subtype and has been well studied in breast cancer (Harari and Yarden, 2000). Hence our Gecco+ not only outperforms existing methods but also selects genes which are relevant to biology and have been implicated in previous scientific studies.
Next we evaluate the performance of iGecco+ for mixed multi-view data sets and investigate the features selected by iGecco+.

Multi-omics Data
One promising application for integrative clustering for multi-view data lies in integrative cancer genomics. Biologists seek to integrate data from multiple platforms of highthroughput genomic data to gain a more thorough understanding of disease mechanisms and detect cancer subtypes. In this case study, we seek to integrate four different types of genomic data to study how epigenetics and short RNAs influence the gene regulatory system in breast cancer.
We use the data set from The Cancer Genome Atlas Network (2012). Lock and Dunson (2013) analyzed this data set using integrative methods and we follow the same data pre-processing procedure: i) filter out genes in expression data whose standard deviation is less than 1.5, ii) take square root of methylation data, and iii) take log of miRNA data. We end up with a data set of 348 tumor samples including: • RNAseq gene expression (GE) data for 645 genes, • DNA methylation (ME) data for 574 probes, • miRNA expression (miRNA) data for 423 miRNAs, • Reverse phase protein array (RPPA) data for 171 proteins.
The data set contains samples used on each platform with associated subtype calls from each technology platform as well as integrated cluster labels from biologists. We use the integrated labels from biologists as true label to compare different methods. Also we merged the subtypes 3 and 4 in the integrated labels as those two subtypes correspond to Luminal A and Luminal B respectively from the predicted label given by gene expression data (PAM50 mRNA). Figure 6 in Appendix H gives the distribution of data from different platforms. For our iGecco+ methods, we use Euclidean distances for gene expression data and protein data as the distributions of these two data sets appear Gaussian; binomial deviances for methylation data as the value is between [0, 1]; Manhattan distances for miRNA data as the data is highly-skewed.
We compare our adaptive iGecco+ with other existing methods. From Table 10, we see that our method outperforms all the existing methods.
We also investigate the features selected by adaptive iGecco+, shown in Table 11, and find that our method is validated as most are known in the breast cancer literature. For example, FOXA1 is known to segregate the luminal subtypes from the others (Badve et al., 2007), and AGR3 is a known biomarker for breast cancer prognosis and early breast cancer detection from blood (Garczyk et al., 2015). Several well-known miRNAs were selected including MIR-135b, which is upregulated in breast cancer and promotes cell growth (Hua et al., 2016) as well as MIR-190 which suppresses breast cancer metastasis (Yu et al., 2018). Several known proteins were also selected including ERalpha, which is overexpressed in early stages of breast cancer (Hayashi et al., 2003) and GATA3 which plays an integral role in breast luminal cell differentiation and breast cancer progression (Cimino-Mathews et al., 2013).
We also visualize the resulting clusters from adaptive iGecco+. Figure 4 shows the cluster heatmap of multi-omics TCGA data with row orders determined by cluster assignments from iGecco+ and left bar corresponding to the integrated cluster labels from biologists. We see that there is a clear separation between groups and adaptive iGecco+ identifies meaningful subtypes. The black bars at the bottom of each data view correspond to the selected features in Table 11.

Discussion
In this paper, we develop a convex formulation of integrative clustering for high-dimensional mixed multi-view data. We propose a unified, elegant methodological solution to two critical issues for clustering and data integration: i) dealing with mixed types of data and ii) selecting sparse, interpretable features in high-dimensional settings. Specifically, we show that clustering for mixed, multi-view data can be achieved using different data-specific convex losses with a joint fusion penalty. We also introduce a shifted group-lasso penalty that shrinks noise features to their loss-specific centers, hence selecting features that play important roles in separating groups. In addition, we make an optimization contribution by proposing and proving the convergence of a new general multi-block ADMM algorithm with sub-problem approximations that efficiently solves our problem. Empirical studies show that iGecco+ outperforms existing clustering methods and selects sparse, interpretable features in separating clusters.
This paper focuses on the methodological development for integrative clustering and feature selection, but there are many possible avenues for future research related to this work. For example, we expect in future work to be able to show that our methods inherit the strong statistical and theoretical properties of other convex clustering approaches such as clustering consistency and prediction consistency. An important problem in practice is choosing which loss function is appropriate for a given data set. While this is beyond the scope of this paper, an interesting direction for future research would be to learn the appropriate convex loss function in a data-driven manner. Additionally, many have shown block missing structure is common in mixed data (Yu et al., 2019;Xiang et al., 2013). A potentially interesting direction for future work would be to develop an extension of iGecco+ that can appropriately handle block-missing multi-view data. Moreover, Weylandt et al. (2020) developed a fast algorithm to compute the entire convex clustering solution path and used this to visualize the results via a dendrogram and pathwise plot. In future work, we expect that algorithmic regularization path approaches can also be applied to our method to be able to represent our solution as a dendrogram and employ other dynamic visualizations. Finally, while we develop an efficient multi-block ADMM algorithm, there may be further room to speed up computation of iGecco+, potentially by using distributed optimization approaches.
In this paper, we demonstrate that our method can be applied to integrative genomics, yet it can be applied to other fields such as multi-modal imaging, national security, online advertising, and environmental studies where practitioners aim to find meaningful clusters and features at the same time. In conclusion, we introduce a principled, unified approach to a challenging problem that demonstrates strong empirical performance and opens many directions for future research.
Proposition 1 and 2 are direct extensions from Proposition 2.1 of Chi and Lange (2015). Notice they proved the solution path depends continuously on the tuning parameter γ and the weight matrix w. It follows that the argument can be also applied to tuning parameter α, the loss weight π k , and feature weight ζ j k . Also it is obvious that the loss ℓ ⋅ is continuous with respect to the data, X.
We show in detail how to prove Proposition 3 in the following. First we rewrite the objective F γ, α U as: By definition, loss-specific cluster center is x k = argmin u i 1 n ℓ k (X i k u). Since ℓ k is convex, it is equivalent to u such that ∂ i ℓ k (X i k u) = 0. Hence, ∂ i ℓ k (X i k x k ) = 0.
We use the similar proof approach of Chi and Lange (2015). A point X furnishes a global minimum of the convex function F X if and only if all forward directional derivatives d Θ F X at X are nonnegative. Here Θ = Θ 1 , ⋯, Θ K represents a direction in the space of possible concatenated centroids, where Θ k ∈ ℝ n × p k . We calculate the directional derivatives: Hence, we have: Take γ sufficiently large such that: When all w ii′ > 0, one can take any γ that exceeds K ⋅ 2 n max i, i′, k π k ‖ ∂ℓ k (X i . k , x k )‖ 2 w ii′ .
For any pair i and i′ there exists a path i k ⋯ l i′ along which the weights are positive. It follows that We have Hence the forward directional derivative test is satisfied for any γ ≥ n 2 β.
On the other hand, for fixed γ, the generalized Cauchy-Schwartz inequality implies: Hence, we have: Take α sufficiently large so that When all ζ j k > 0, one can take any α that exceeds max j, k π k ‖ ∂ℓ k (X . j k , x j k ⋅ 1 n ‖ 2 ζ j k . In general, set α ≥ 1 min ζ j k > 0 ζ j k ⋅ max j, k π k ‖ ∂ℓ k (X . j k , x j k ⋅ 1 n )‖ 2 . It is easy to check the forward directional derivative test is satisfied. ■

Appendix B.: Proof of Theorem 4
Recall the iGecco+ problem is: We can recast the original iGecco+ problem as the equivalent constrained problem: where ℓ k refers to the differentiable losses and ℓ k′ refers to the non-differentiable losses. We use multi-block ADMM algorithm (Algorithm 6) to solve the problem above.
To prove convergence of Algorithm 5, we first show that multi-block ADMM Algorithm 6 converges to a global minimum. Then we show that we can proximal-linearize the sub-problems in the primal updates of Algorithm 6 with proved convergence and this is equivalent to Algorithm 5. Without loss of generality, we assume we have one differentiable loss ℓ 1 ⋅ and one non-differentiable distance-based loss ℓ 2 ⋅ .

Algorithm 6
Multi-block ADMM to solve iGecco+ while not converged do for all k = 1, ⋯, K do U k = argmin U π k ℓ k (X k , U) + ρ 2 ‖DU − V k + Λ k ‖ F 2 + α j 1 p k ζ j k ‖U j x j k 1 n ‖ 2 U (k′) = argmin To prove convergence of Algorithm 6, we first propose a sufficient condition for the convergence of four-block ADMM and prove it holds true. This is an extension of the convergence results in Section 2 of the work by Chen et al. (2016). Suppose the convex optimization problem with linear constraints we want to minimize is min θ 1 x 1 + θ 2 x 2 + θ 3 x 3 + θ 4 x 4 s.t. A 1 x 1 + A 2 x 2 + A 3 x 3 + A 4 x 4 = b .
The multi-block ADMM has the following form. Note here, the superscript in x i k + 1 refers to the (k + 1) th iteration in the ADMM updates. We have: where We establish Lemma 5, a sufficient condition for convergence of four-block ADMM: Lemma 5 (Sufficient Condition for Convergence of Four-block ADMM) A sufficient condition ensuring the convergence of (4) to a global solution of (3) is: A 2 T A 3 = 0, A 2 T A 4 = 0, A 3 T A 4 = 0.
We prove Lemma 5 at the end of this section. Note that Lemma 5 is stated in vector form and therefore we need to transform the constraints in the original iGecco+ problem (2) from matrix form to vector form in order to apply Lemma 5. Note that )vec(U k T ) = vec(V k T ). Hence we can write the constraints in (2) as: Next we show that the constraints in (2)  According to the first-order optimality conditions of the minimization problems in (4), we have: which is also the first-order optimality condition of the scheme: } (x 2 k + 1 , x 3 k + 1 , x 4 k + 1 ) = argmin{θ 2 x 2 + θ 3 x 3 + θ 4 x 4 − (λ k ) T A 2 x 2 + A 3 x 3 + A 4 x 4 Clearly, (5) is a specific application of the original two-block ADMM to (3) by regarding (x 2 , x 3 , x 4 ) as one variable, [A 2 , A 3 , A 4 ] as one matrix and θ 2 (x 2 ) + θ 3 (x 3 ) + θ 4 (x 4 ) as one function. ■

Appendix C.: Gecco+ for Differentiable Losses
In this section, we propose algorithms to solve Gecco+ when the loss ℓ is differentiable and gradient is Lipschitz continuous. In this case, we develop a fast two-block ADMM algorithm without fully solving the U sub-problem. Our result is closely related to the proximal linearized ADMM literature Lu et al., 2016). Also solving the sub-problem approximately is closely connected with the generalized ADMM literature (Deng and Yin, 2016).
In the following sections, we discuss algorithms to solve Gecco+ instead of iGecco+ for notation purposes as we would like to include iteration counter indices in the algorithm for illustrating backtracking; but we can easily extend the algorithm to solve iGecco+. To begin with, we clarify different notations in Gecco+ and iGecco+: the superscript in U (k) in iGecco+ refers to the k th data view while U (k) in Gecco+ refers to the k th iteration counter in the ADMM updates. We omit iteration counter indices in all iGecco+ algorithm for notation purposes and use the most recent values of the updates.

C.1 Two-block ADMM in Matrix Form
minimize U, V ℓ X, U + γ l ε w l ‖V l ‖ 2 P 1 V; w + α j 1 p ζ j ‖U j x j 1 n ‖ 2 subject to DU − V = 0 .
Like in convex clustering (Chi and Lange, 2015;Weylandt et al., 2020), we index a centroid pair by l = l 1 , l 2 with l 1 < l 2 , define the set of edges over the non-zero weights ε = l = l 1 , l 2 : w l > 0 , and introduce a new variable V l . = U l 1. − U l 2. to account for the difference between the two centroids. Hence V is a matrix containing the pairwise differences between connected rows of U. Also D is the directed difference matrix corresponding to the non-zero fusion weights defined in the work of Weylandt et al. (2020).
We can show that the augmented Lagrangian in scaled form is equal to: where the dual variable is denoted by Λ.
To update U, we need to solve the following sub-problem: Let U = U − X. The sub-problem becomes: where u j is the j th column of U. For each ADMM iterate, we have: . This can be solved by running iterative proximal gradient to full convergence: Here U (k,m) refers to the m th inner iteration counter in the U sub-problem out of the k th outer iteration counter of the ADMM update. It is straightforward that this is computationally expensive. To address this, we propose to solve the U sub-problem approximately using just a one-step proximal gradient update and prove convergence in the next section. This approach is based on proximal linearized ADMM Lu et al., 2016), which solves the sub-problems efficiently by linearizing the differentiable part and then applying proximal gradient due to the non-differentiable part. To ensure convergence, the algorithm requires that gradient should be Lipschitz continuous. The V and Λ updates are just the same as in regular convex clustering.
We adopt such an approach and develop the proximal linearized 2-block ADMM (Algorithm 7) to solve Gecco+ when the loss is differentiable and gradient is Lipschitz continuous.

Algorithm 7
Proximal linearized 2-block ADMM when the loss is differentiable and gradient is Lipschitz continuous -matrix form Further, if the U sub-problem can be decomposed to p separate U .j sub-problems where the augmented Lagrangian for each now is a sum of a differentiable loss, a quadratic term and a sparse group-lasso penalty, we propose to use proximal gradient descent for each separate U .j sub-problem. In this way, we yield adaptive step size for each U .j sub-problem and hence our algorithm enjoys better convergence property than updating U's all together. (In the latter case, the step size becomes fairly small as we are moving all U to some magnitude in the direction of negative gradient.) To achieve this, we assume that the loss is elementwise, which means we can write the loss function as a sum of p terms. (The loss can be written as i ℓ X i U i j ℓ X . j , U . j = i j q x ij u ij where q(·) is the element-wise version of the loss while ℓ(·) is the vector-wise version of the loss.) We see that every deviance-based loss satisfies this assumption. Moreover, by decomposing to p sub-problems, we can solve

C.2 Two-block ADMM in Vector Form in Parallel
Suppose the U sub-problem can be decomposed to p separate U .j sub-problems mentioned above. The augmented Lagrangian now becomes: In this way we can perform block-wise minimization. Now minimizing the augmented Lagrangian over U is equivalent to minimizing over each U . j , j = 1, ⋯, p: minimize U . j ℓ(X . j , U . j ) + ρ 2 ‖DU . j − V . j + Λ . j ‖ 2 2 + αζ j ‖U . j − x j ⋅ 1 n ‖ 2 .
Let u j = U . j − x j ⋅ 1 n . The problem above becomes: Similarly, this can be solved by running iterative proximal gradient to full convergence. However, as mentioned above, we propose to solve the U sub-problem approximately by taking just a one-step proximal gradient update and prove convergence. Still this approach is based on proximal linearized ADMM Lu et al., 2016), which solves the sub-problems efficiently by linearizing the differentiable part and then applying proximal gradient due to the non-differentiable part. To ensure convergence, the algorithm requires that gradient should be Lipschitz continuous.
We propose Algorithm 8 to solve Gecco+ when ℓ is differentiable and gradient is Lipschitz continuous in vector form. Note the U-update in Algorithm 8 is a just a vectorized version of that in Algorithm 7 if we use fixed step size s k for each feature j. We use the vector form update here since it enjoys better convergence property mentioned above and we use this form to prove convergence. Next we prove the convergence of Algorithm 8.

C.3 Proof of Convergence
Theorem 6 If ℓ is convex and differrentiable and ∇ℓ is Lipschitz continuous, then Algorithm 8 converges to a global solution. Further, if ℓ is strictly convex, it converges to the unique global solution.
Proof: We will show that the U-update in Algorithm 8 is equivalent to linearizing the U sub-problem and then applying a proximal operator, which is proximal linearized ADMM.
Note that each U .j sub-problem is: For simplicity of notation, we replace U .j with u j in the following: u j k + 1 = argmin u j ℓ(X . j , u j ) + ρ 2 ‖Du j − V . j k + Λ . j k ‖ 2 2 + αζ j ‖u j − x j ⋅ 1 n ‖ 2 .
Recall the definition of proximal operator: Therefore, the u j update is just a proximal gradient descent update: u j k + 1 = prox s k ⋅ αζ j ‖ ⋅ ‖ 2 (u j k − x j ⋅ 1 n − s k ⋅ [ ∇ℓ(X . j , u j k ) + ρD T (Du j k − V . j k + Λ . j k )]) . Now we plug back and get the u j , i.e., U . j update: which is equivalent to the U .j update in Algorithm 8.
The V and Λ update is just the same as the one in convex clustering. Hence Algorithm 8 is equivalent to the form of proximal linearized ADMM by Liu et al. (2013);Lu et al. (2016) and hence converges to a global solution as long as ∇ℓ is Lipschitz continuous. Note that Algorithm 7 is equivalent to Algorithm 8 with a fixed step size s k . (We can choose s k to be the minimum step size s k over all features.) Therefore, Algorithm 7 also converges to a global solution as long as ∇ℓ is Lipschitz continuous. Further, if ℓ is strictly convex, the optimization problem has unique minimum and hence Algorithm 7 and 8 converges to the global solution. ■ Note: • In Liu et al. (2013); Lu et al. (2016), the algorithm requires that ∇ℓ is Lipschitz continuous to guarantee convergence. We know that ∇ℓ being Lipschitz continuous is equivalent to ℓ being strongly smooth. It is easy to show that the Hessian of log-likelihood of exponential family and GLM deviance is upper bounded since in (generalized) convex clustering, the value of U is bounded as it moves along the regularization path from X to the loss-specific center; also to avoid numerical issues, we add trivial constraint that u ij > 0 when zero is not defined in the log-likelihood/deviance. Hence the condition for convergence of proximal linearized ADMM is satisfied.

•
To obtain a reasonable step size s k , we need to compute the Lipschitz constant.
However, it is non-trivial to calculate the Lipschitz constant for most of our general losses. Instead, we suggest using backtracking line search procedure proposed by Beck and Teboulle (2009); Parikh et al. (2014), which is a common way to determine step size with guaranteed convergence in optimization.
Empirical studies show that choosing step size with backtracking in our framework also ensures convergence. The details for backtracking procedure are discussed below.

C.4 Backtracking Criterion
In this section we discuss how to choose the step size s k in Algorithm 8. As mentioned, usually we employ a fixed step size by computing the Lipschitz constant as in the squared error loss case; but in our method, it is hard to compute the Lipschitz constant for most of our general losses. Instead, we propose using backtracking line search procedure proposed by Beck and Teboulle (2009); Parikh et al. (2014), which is a common way to determine step size with guaranteed convergence in optimization.
We adopt the backtracking line search procedure proposed by Beck and Teboulle (2009);Parikh et al. (2014). At each iteration, while g(u j − tG t (u j )) > g(u j ) − t ∇g(u j ) T G t (u j ) + t 2 ‖G t (u j )‖ 2 2 i.e., g(prox t (u j − t ∇g(u j ))) > g(u j ) − ∇g(u j ) T (u j − prox t (u j − t ∇g(u j ))) We still adopt the one-step approximation and hence suggest taking a one-step proximal update with backtracking to solve the U sub-problem. To summarize, we propose Algorithm 9, which uses proximal linearized 2-block ADMM with backtracking when the loss is differentiable and gradient is Lipschitz continuous.

C.5 Alternative Algorithm for Differentiable Losses
It should be pointed out that there are many other methods to solve the U sub-problem when the loss ℓ is differentiable. We choose to use proximal gradient descent algorithm as there is existing literature on approximately solving the sub-problem using proximal gradient under ADMM with proved convergence Lu et al., 2016). But there are many other optimization techniques to solve the U sub-problem such as ADMM.
In this subsection, we show how to apply ADMM to solve the U sub-problem and specify under which conditions this method is more favorable. Recall to update U, we need to solve the following sub-problem: We use ADMM to solve this minimization problem and can now recast the problem above as the equivalent constrained problem: minimize U, V, Λ, R j 1 p ℓ(X j U j ) + ρ 2 ‖DU − V + Λ‖ F 2 + α j 1 p ζ j ‖r j ‖ 2 P 2 R; ζ subject to U−X = R .

Wang and Allen Page 38
Proximal linearized 2-block ADMM with backtracking when the loss is differentiable and gradient is Lipschitz continuous Input: X, γ, w, α, ζ The augmented Lagrangian in scaled form is: where the dual variable for V is denoted by Λ; the dual variable for R is denoted by N.
The U .j sub-problem in the inner nested ADMM is: Wang and Allen Page 39 Now, we are minimizing a sum of a differentiable loss ℓ and two quadratic terms which are all smooth. Still the U .j sub-problem does not have a closed-form solution for general convex losses and we need to run an iterative descent algorithm (such as gradient descent, Newton method) to full convergence to solve the problem. Similarly, to reduce computation cost, we take a one-step update by applying linearized ADMM (Lin et al., 2011) to the U .j

Algorithm 10
Full-step multi-block algorithm for Gecco+ with Bernoulli log-likelihood ℓ k : end while V (k) = prox γ/ρP 1 ( ⋅ ; w) (DU (k) + Λ (k − 1) ) Λ (k) = Λ (k − 1) + (DU (k) − V (k) ) end while sub-problem. The U .j update in the inner ADMM now becomes: In this case, empirical studies show that taking a one-step Newton update is favored than a one-step gradient descent update as the former enjoys better convergence properties and generally avoids backtracking. However, inverting a Hessian matrix is computationally burdensome at each iteration when n is large. Exceptions are for Euclidean distances case where there is a closed-form solution for the U .j update and for Bernoulli log-likelihood case where the Hessian of the loss can be upper bounded by a fixed matrix. In the latter case, we propose to pre-compute the inverse of that fixed matrix instead of inverting a Hessian matrix at each iteration. To illustrate this, we write out the U .j sub-problem of Gecco+ with Bernoulli log-likelihood: U . j = argmin U . j i 1 n x ij u ij log 1 e u ij + ρ 2 ‖DU . j − V . j + Λ . j ‖ F 2 The Hessian is diag e u ij 1 + e u ij 2 + ρD T D + ρI which can be upper bounded by 1 4 I + ρD T D + ρI.
We propose to replace the Hessian with this fixed matrix in Newton method and use its inverse. This is closely related to the approximate Hessian literature (Krishnapuram et al., 2005;Simon et al., 2013). In this way, we just pre-compute this inverse matrix instead of inverting the Hessian matrix at each iteration, which dramatically saves computation. We give Algorithm 10 to solve Gecco+ for Bernoulli log-likelihood with a one-step update to solve the U sub-problem. Empirical studies show that this is faster than taking the inner nested proximal gradient approach as we generally don't need to perform the backtracking step.
Yet, Algorithm 10 is slow as we need to run iterative inner nested ADMM updates to full convergence. To address this, as mentioned, we can take a one-step update of the inner nested iterative ADMM algorithm. To see this, we can recast the original Gecco+ problem as: minimize U k , V ℓ X, U + γ l ε w l ‖V l ‖ 2 P 1 V; w + α j 1 p ζ j ‖r j ‖ 2 P 2 R;ζ subject to DU−V = 0, U − X = R .
We apply multi-block ADMM to solve this optimization problem and hence get Algorithm 11. As discussed above, we take a one-step update to solve the U sub-problem with linearized ADMM and use the inverse of fixed approximate Hessian matrix, M 1 .

Appendix D.: Gecco+ for Non-Differentiable Losses
In this section, we propose an algorithm to solve Gecco+ when the loss ℓ is nondifferentiable. In this case, we develop a multi-block ADMM algorithm to solve Gecco+ and prove its algorithmic convergence.

D.1 Gecco+ Algorithm for Non-Differentiable Losses
f Z = j 1 p ‖z j ‖ 1 ‖vec Z ‖ 1 and g X, U = X − U. We specify the affine function g as we want to augment the non-differentiable term in the loss function ℓ.
We can rewrite the problem as: We can now recast the problem above as the equivalent constrained problem: where X is an n × p matrix with the j th column equal to scalar x j .
The augmented Lagrangian in scaled form is: where the dual variable for V is denoted by Λ; the dual variable for Z is denoted by Ψ; the dual variable for R is denoted by N.
Since we assume g to be affine, i.e., g(X, U) = AX + BU + C, the augmented Lagrangian in scaled form can be written as: It can be shown that the U sub-problem has a closed-form solution.
Note that, hinge loss is also non-differentiable and we can write g(X, U) = 1 − U ο X where For distance-based losses, the loss function can always be written as: ℓ(X, U) = f(X − U), which means g(X, U) = X − U. Then the augmented Lagrangian in scaled form can be simplified as: Now the U sub-problem has a closed-form solution: This gives us Algorithm 13 to solve Gecco+ for non-differentiable distance-based loss.

Algorithm 13
Multi-block ADMM for non-differentiable distance-based loss while not converged do Algorithm 13 can be used to solve Gecco+ problem with various distances such as Manhattan, Minkowski and Chebychev distances by applying the corresponding proximal operator in the Z-update. For example, for Gecco+ with Manhattan distances, the Z-update is just applying element-wise soft-thresholding operator. For Gecco+ with Chebychev distances, the proximal operator in the Z-update can be computed separately across the rows of its argument and reduces to applying row-wise proximal operator of the infinity-norm. For Gecco+ with Minkowski distances, we similarly apply row-wise proximal operator of the ℓ q -norm.
Next we prove the convergence of Algorithm 13. Wang and Allen Page 44 Like before, we can rewrite the problem as: We can now recast the problem above as the equivalent constrained problem: Here, f(Z) = max(0, Z). With a slight abuse of notation, we refer f to applying element-wise maximum to all entries in the matrix. We set g(X, U) = 1 − U ο X where 1 is a matrix of all one and "ο" is the Hadamard product. X is an n × p matrix with the j th column equal to scalar x j .
The augmented Lagrangian in scaled form is: The U sub-problem now becomes: The first-order optimality condition is: which can be written as: To solve U from the above equation, one way is to first find the SVD of the leading coefficient: From this decomposition we create two sets of diagonal matrices: The Hadamard product can now be replaced by a sum Here we have to compute the pseudo-inverse of a matrix which is computationally expensive in practice. To avoid this, we adopt generalized ADMM approach proposed by Deng and Yin (2016) where the U sub-problem is augmented by a positive semi-definite quadratic operator. In our case, our modified U sub-problem becomes: The first-order optimality condition now becomes: We have: where H = (D T D + I + 1). Hence we have analytical update: It is easy to see that the V, Z and R updates all have closed-form solutions.

Appendix E.: Multinomial Gecco+
In this section, we briefly demonstrate how Gecco+ with multinomial losses is formulated, which is slightly different from the original Gecco+ problems. Suppose we observe categorical data as follows (K = 3): X n × p = 1 1 1 2 2 2 3 3 3 .
We can get the indicator matrix X (k) for each class k as: Then we concatenate X (1) , X (2) , X (3) and get X n × p * K = X 1 X 2 X 3 . This is equivalent to the dummy coding of the categorical matrix X after some row/column shuffle: It is obvious that measuring the difference of two observations by comparing rows of X is better than simply comparing the Euclidean distances of rows of original data matrix X. Also parameterizing in X is beneficial for computing the multinomial log-likelihood or deviance. Hence we concatenate all columns of X (k) as input data. Similarly, we concatenate all columns of the corresponding U (k) and then fuse U in a row-wise way.

E.1 Gecco with Multinomial Log-likelihood
Gecco with multinomial log-likelihood can be formulated as: Wang and Allen Page 48 where x ijk refers to the elements of indicator matrix X ij k discussed previously and

E.2 Gecco+ with Multinomial Log-likelihood
Gecco+ with multinomial log-likelihood can be formulated as: ⋮ u njk and x j k is the loss-specific center for j variable in k th class.

Appendix F.: Loss-specific Center Calculation
In this section, we show how to calculate the loss-specific center in Table 1.

F.1 Continuous Data
For continuous data, we consider Gecco with Euclidean distances.

F.2 Count Data
For count-valued data, we consider Gecco with Poisson log-likelihood/deviance, negative binomial log-likelihood/deviance and Manhattan distances. Taking derivative, we get: i 1 n x ij + nexp(u j ) = 0 exp(u j ) = x j u j = log(x j ) u = log(x) .

F.2.2 POISSON DEVIANCE
Let u be the fusion vector and u = (u 1 , ⋯ , u p ). The problem above becomes: Taking derivative, we get:  When total fusion occurs, u ij = u i′j , ∀i ≠ i′. Let u be the fusion vector and u = (u 1 , ⋯, u p ). The problem above becomes: Taking derivative, we get:

F.2.4 NEGATIVE BINOMIAL DEVIANCE
The formulation above is equivalent to: Let u be the fusion vector and u = u 1 , ⋯, u p . The problem above becomes:

F.2.5 MANHATTAN DISTANCE
When total fusion occurs, U i . = U i′ . , ∀i ≠ i′. Let U i . = U i′ . = u. We have: For each j, we have: We know that the optimal u j is just the median of x ij for each j.

F.3 Binary Data
For binary data, we consider Gecco with Bernoulli log-likelihood, binomial deviance and hinge loss.

F.3.1 BERNOULLI LOG-LIKELIHOOD
When total fusion occurs, u ij = u i′j , ∀i ≠ i′. Let u be the fusion vector and u = u 1 , ⋯, u p . The problem above becomes: minimize U i n j p x ij u j + log(1 + e u i ) .

F.3.2 BINOMIAL DEVIANCE
Let u be the fusion vector and u = u 1 , ⋯, u p . The problem above becomes: Taking derivative, we get: Let u be the fusion vector and u = u 1 , ⋯, u p . The problem above becomes:

F.4 Categorical Data
For categorical data, we consider Gecco with multinomial log-likelihood and deviance.

F.4.1 MULTINOMIAL LOG-LIKELIHOOD
When total fusion occurs, u ijk = u i′jk , ∀i ≠ i′. Let u be the fusion vector and u k = u 1k , ⋯, u pk . The problem above becomes: Taking derivative with respect to u jk , we get: When total fusion occurs, u ijk = u i′jk , ∀i ≠ i′. Let u be the fusion vector and u k = u 1k , ⋯, u pk . The problem above becomes: We can write the constraint in Lagrangian form: Taking derivative with respect to u jk , we get: We have: Therefore, u jk = i 1

Figure 6:
Histograms of data from different platforms for multi-omics TCGA data set. Both gene expression data and protein data appear Gaussian; methylation data is proportion-valued; miRNA data is highly-skewed.

Appendix I.: Null Deviance
In this appendix, we list the null deviance D for some common losses used in iGecco. Recall in the iGecco formulation, π k , which are set inversely proportional to the null deviance evaluated at the loss-specific center, are scaling factors to ensure that the losses are measured at the same scale in the objective function.
By definition, the null deviance D evaluated at the loss-specific center, refers to ℓ(X, X) where each j th column of X is the loss-specific center x j k for that column/feature. Table 12 shows the null deviance D for some common losses.  Here x j refers to the mean for the j th column/feature while median(x j ) refers to the median for the j th column. Simple calculation shows that the Bernoulli log-likelihood is equivalent to binomial deviance by noting logit x j = log

Appendix J.: Choice of Tuning Parameters
In this appendix, we propose two different approaches to select tuning parameters γ and α in the iGecco+ problem; γ controls the number of clusters while α controls the number of features selected. We first consider stability selection based approach, which has been shown to enjoy nice statistical properties such as selection consistency (Wang, 2010;Fang and Wang, 2012) and been adopted in practice with strong empirical performance . Next, to reduce computation, we consider information criterion based approaches. Finally, we demonstrate empirical results when the number of clusters and features are not fixed but estimated based on the data.

J.1 Stability Selection Based Approach
We first adopt the stability selection based approach for tuning parameter selection and follow closely the approach described in the work of Wang (2010); Fang and . We choose the stability selection based approach because i) its selection consistency has been established and ii) Wang et al. (2018) adopted similar approach for tuning parameter selection and demonstrated strong empirical performance. The rationale behind stability selection is that a good clustering algorithm with optimal tuning parameter should produce clusterings that do not vary much with respect to a small perturbation to the training samples.
By Wang (2010); Fang and Wang (2012), a clustering ψ(x) is defined as a mapping ψ : ℝ p {1, …, q} where q is the given number of clusters. Here, we use q as the number of clusters since we have already used k to represent the k th data view X (k) . A clustering algorithm Ψ(·; q) with a given number of clusters q ≥ 2 yields a clustering mapping ψ(x) when applied to a sample. Here, we choose the number of clusters q, which is equivalent to choosing optimal γ since we can yield the cluster assignment and corresponding q for each γ. To account for the tuning parameter α for feature selection, we further denote the clustering algorithm as Ψ(·; q, α). Then, in our case, for any given pair of q and α, two clustering results can be obtained from two sets of bootstrapped samples. Then the clustering distance d, defined by Fang and , can be computed to measure the dissimilarity between two clustering results. We repeat the procedure multiple times and the optimal tuning parameter pair (q, α) is the one that minimizes the average clustering instability. This gives Algorithm 14, which uses stability selection to choose tuning parameters.

Algorithm 14
Choice of number of clusters q and feature penalty γ in Gecco+/iGecco+ using stability 4. Finally, the optimal number of clusters q and feature penalty α can be estimated by q, α = argmin 2 ≤ q ≤ Q, α S B Ψ, q, α .

J.2 Information Criterion Based Approach
While choosing number of clusters based on stability selection has been shown to enjoy nice properties such as selection consistency (Wang, 2010), such methods are always computationally burdensome. To address this, on the other hand, information criterion based approach has been proposed for tuning parameter selection in convex clustering (Tan and Witten, 2015). We adopt the similar approach and propose the Bayesian information criterion (BIC) based approach to choose optimal number of clusters and features; this gives Algorithm 15.
The criterion to minimize is inspired by the BIC approach for convex clustering proposed by Tan and Witten (2015). Notice that we use different criterion in step 1 and 2 to choose number of clusters and features respectively. When choosing the number of clusters, we use the same loss function as in iGecco since all the features are selected. When choosing the number of features, for the presence of noise features, we find that the weighted loss function with respect to each feature in step 2 works better, as the criterion in step 1 downweights the contribution of signal features by adding informative terms along with the noise terms in the denominator. (When a feature is not selected, ℓ k (X . j k , U . j k ) equals ℓ k (X . j k , X . j k ); when a feature is selected, ℓ k (X . j k , U . j k ) is less than ℓ k (X . j k , X . j k ); if we use the unweighted criterion in step 1, the noise terms will dominate this criterion given a large number of noise features.) The degrees of freedom for choosing number of clusters in step 1 is q while in step 2, the degrees of freedom is sq since we need to estimate q number of cluster centroids for each selected feature.
Stability selection is known to be more stable in terms of choosing number of clusters, with selection consistency theoretically established in literature. Meanwhile, BIC works better in choosing number of features in practice and saves much more computation compared with stability selection. Hence, to take full advantage of both approaches, we propose a sequential tuning parameter selection procedure, outlined in Algorithm 16. We demonstrate the clustering and variable selection accuracy results in Appendix J.3 using the proposed two tuning parameter selection approaches. 3. Update the feature weights ζ j k = 1/‖U . j k − x j k ⋅ 1 n ‖ 2 and fusion weights w ij = I ij k exp(−ϕd(X i . , X i′ . )), where d(X i . , X i′ . ) = k 1 K j 1 p k ‖U j k x j k 1 n ‖ 2 max j ‖U j k x j k 1 n ‖ 2 Fit adaptive iGecco+ with ζ , w and a sequence of γ and α; Find optimal γ and α using the tuning parameter selection procedure in Algorithm 15 or 16.

J.3 Empirical Results in Simulation and Real Data Example
In this section, we demonstrate the empirical results when the number of clusters and features are not fixed but estimated based on the data. Specifically, we apply the two tuning parameter selection schemes proposed above: i) BIC based approach (Algorithm 17 + 15) and ii) stability selection + BIC based approach (Algorithm 17 + 16). We show the results of the tables in Section 4 and 5 when the number of clusters and features are estimated and not fixed to be the oracle. For simulation studies, we show the results of iGecco and iGecco+; Gecco and Gecco+ are special cases of these two. Moreover, we apply the proposed two tuning parameter selection approaches to the three real data examples in Section 5.
For iGecco, the stability selection based approach simplifies to choosing the number of clusters q in Algorithm 14; the BIC based approach refers to using the criterion in step 1 of Algorithm 15 to choose the number of clusters.
For iGecco+, we show the overall F1-score and number of selected features across all three data views. Recall that each data view has 10 true features and hence there are 30 true features in total. Also, the optimal number of clusters is q = 3. We compare the results with iClusterPlus when the number of clusters and features are not fixed. Note we only include estimated number of clusters for iClusterPlus as there is no tuning parameter for the number of selected features. (The authors mentioned feature selection could be achieved by selecting the top features based on Lasso coefficient estimates.) Table 13 and 14 show the results of Table 3, 4 and 5 in Section 4.2 when the number of clusters (and features in iGecco+) is not fixed but estimated based on the data. Table 13 and 14 suggest that our proposed BIC based approach selects the correct number of clusters and features most of the time. On the other hand, information criterion based approaches save much more computation than stability selection. Hence, we recommend the BIC based approach for choosing tuning parameters which demonstrates strong empirical performance and saves computation. Yet, one might have their own justified choice of approach or information criterion to select tuning parameters. Adjusted Rand index and estimated number of clusters of iGecco and iCluster for mixed multi-view data of Table 3 in Section 4.2 when the number of clusters is not fixed but estimated based on the data.  Adjusted Rand index, F 1 score, along with estimated number of clusters and features of adaptive iGecco+ and iClusterPlus for high-dimensional mixed multi-view data of Table 4 and 5 in Section 4.2 when the number of clusters and features are not fixed but estimated based on the data. We only include estimated number of clusters for iClusterPlus as there is no tuning parameter for the number of selected features. Also, we apply the proposed two tuning parameter selection approaches to the real data examples in Section 5. Table 15 shows the estimated number of clusters. Again, the BIC based approach selects the correct number of clusters for all three cases. Adjusted Rand index, along with estimated number of clusters using adaptive (i)Gecco+ for real data of Table 6, 8 and 10 in Section 5 when the number of clusters and features are not fixed but estimated based on the data. For the first two data set, we show the results of Manhattan Gecco+.

Appendix K.: Stable to Perturbations of Data
In this appendix, we demonstrate that clustering assignments of iGecco+ are stable to perturbations in the data as shown in Proposition 2 of Section 2.4.
To show this, we include a simulation study similar to the one by Chi et al. (2017). We first apply adaptive iGecco+ on the original data to obtain baseline clustering. Then we add i.i.d. noise to each data view to create a perturbed data set on which we apply the same iGecco+ method. Specifically, for Gaussian data view, we add i.i.d. N(0, σ 2 ) noise where σ = 0.5, 1.0, 1.5; for count data view, we add i.i.d. Poisson noise; for binary data view, we randomly shuffle a small proportion of the entries. We compute the adjusted Rand index between the baseline clustering and the one obtained on the perturbed data. We adopt the same approach for other existing methods. Table 16 shows the average adjusted Rand index of 10 replicates.
For all values of σ, we see that iGecco+ tends to produce the most stable results. Stability and reproducibility of adaptive iGecco+ on simulated data. Adaptive iGecco+ and other existing methods are applied to the simulated data to obtain baseline clusterings. We then perturb the data by adding i.i.d. noise. Specifically, for Gaussian data view, we add i.i.d. N(0, σ 2 ) noise where σ = 0.5, 1.0, 1.5; for count data view, we add i.i.d. Poisson noise; for binary data view, we randomly shuffle a small proportion of the entries. We compute the adjusted Rand index (ARI) between the baseline clustering and the one obtained on the perturbed data.

Appendix L.: Noisy Data Sources
In this appendix, we show the performance of iGecco+ when a purely noisy data view is observed.
Often in real data, not all the data views observed contain clustering information, i.e., one or more data views might be pure noise. Our iGecco+ is able to filter out noisy data views by adaptively shrinking all the (noise) features in these data sets towards the loss-specific centers. We include a simulation study to demonstrate the performance of iGecco+ in the presence of some purely noisy data sets. Similar to the base simulation, each simulated data set consists of n = 120 observations with 3 clusters. Each cluster has an equal number of observations. Only the first data view contains clustering signal with the first 30 features being informative while the rest features being noisy; the rest two data views are pure noise. Table 17 shows that iGecco+ still performs well in the presence of noise data sources by adaptively shrinking noise features.

Appendix M.: Computation Time
In this section, we provide some computation run time results of iGecco(+) with different sample sizes and dimensions. We include results for both run times per iteration of the ADMM algorithm in Table 18 as well as the full training time for tuning parameter selection in Table 19. All timing results are run on a Dell XPS 15 with a 2.4 GHz Intel i5 processor and 8 GB of 2666 MHz DDR4 memory.
For training, we use BIC based approach to select tuning parameters for both iGecco and iGecco+ as it works well in practice and saves computation. It takes more time to train iGecco+ than iGecco as iGecco+ needs to select two tuning parameters (the number of clusters and features). Yet, our BIC based approach selects optimal tuning parameters in a reasonable amount of time. Note, for sample size n = 120 and feature size p 1 = 200, p 2 = 100, p 3 = 50, it takes iClusterPlus hours for tuning parameter selection (using tune.iClusterPlus function in R).   Simulation results of non-Gaussian data: (S1A) We increase number of noise features for spherical data with outliers; (S2) We increase number of noise features for non-spherical data with outliers; (S3) We increase number of noise features for count-valued data; (S1B) We increase noise level for spherical data with outliers; (S1C) We further increase number of noise features for spherical data with outliers in high dimensions. The adaptive Gecco+ outperforms existing methods in high dimensions.